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Abstract.  Finite  element  problems  can  often  naturally  be  divided  into 
subproblems  which  correspond  to  subregions  into  which  the  region  has 
been  partitioned  or  from  which  it  was  originally  assembled.  A  class  of 
iterative  methods  are  discussed  in  which  these  subproblems  are  solved 
by  direct  methods,  while  the  interaction  across  the  curves  or  surfaces 
which  divide  the  region  is  handled  by  a  conjugate  gradient  method.  A 
mathematical  framework  for  this  work  is  provided  by  regularity  theory 
for  elliptic  finite  element  problems  and  by  block  Gaussian  elimination. 
A  full  development  of  the  theory,  which  shows   that   certain  of   these 
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methods   are   optimal,    is   given  for  Lagrangian  finite  element 
approximations  of  second  order  linear  elliptic  problems  in   the   plane. 
Results  from  numerical  experiments  are  also  reported. 

1.    INTRODUCTION  3ijlo'      sr'r  shvloroo      cr 

In  this  papfit^r  whl^h  s>is  %ft:^J expanded  version  of  Bj^rstad  and 
Widlund  [4],  we  deveiopr'and  aftaly^fe'  soma  iterative  methods  for  the 
solution  of  ell  Iptic  probleite^-dh  a  region  regarded  as  a  union  of 
subregions.  We  ar?.,in^partical4r^- interested  in  the  design  of  methods 
for  which  the  interaction  between  the  subproblems,  i.e.  the  discrete 
elliptic  problems  oni  the'  ■subreglbrii'i  is  computed  by  a  conjugate 
gradient  method,  but  the;  subpr6bleto§'^afe  solved  using  a  direct  method. 

The  partition  of  elliptic --pfobl^m^  into  subproblems  is  a  very 
natural  idea  and  is  very  much  ift^a*e?iM  practice.  In  finite  element 
computations,  the  modeling  of --^Sn-e'6tire  structure  can  profitably  be 
organized  by  discretizing  the  giv^'  pd^ial  differential  equations  by 
finite  elements  on  subregions  into  "which  the  region  is  partitioned  or 
from  which  it  has  originally  been  a&sembled.  The  Cholesky  triangular 
factorization  of  the  stiffness  matrices  of  the  subregions  is  then 
computed.  These  tasks  can  obviously  be  assigned  to  different 
engineering  groups,  computer  systems  or  processors  with  coordination 
required  only  at  the  interfaces  between  the  substructures  to  assure  a 
matching  of  the  finite  element  triangularizations.  The  solution  of  the 
entire  problem  is  completed  by  taking  the  interaction  between  the 
substructures  into  account,  see  Przemieniecki  [28,29].  These  ideas  are 
also  often  used  recursively  and  they  are  particularly  attractive  if 
several   substructures   are   identical   or   if   some  of  them  have  been 
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analyzed  previously.   Such  a  situation  arises   for  example  if   an 
analysis   is  repeated  after  the  redesign  or  damage  of  one  or  relatively 
few  of  the  substructures,  see  Bell  et  aX*     rii]      and  Hatlestad  and 
Mellingen  [18].  ' 

It  is  normal  practice  to  conclude  the  solution  process  by 
computing  the  matrix  which  de^jcxibj^  Slie  niateraptlpn  between  the 
substructures  and  its  triangul^^f^  fag^psrl^atlonor  and-  by  using  the 
triangular  factors  to  solve  the  iineaife  I  systems  :  There  are  many 
interesting  sparse  matrix  and  sof6>/9re-:risa9esi involved  which  are  being 
studied  actively,  see  e.g.  Duff  and  Reidv{i3]'.  Weralso  note  that  the 
important  work  by  George  [14]  and  oth^rsfeon. nested  dissection  of  finite 
element  and  related  systems  ^qf  sit^i^fikE-i  daigebraic  equations  can  be 
regarded  as  a  recursive  applipatiQii^rof  ;Si-Biilar  ideas.  The  results  of 
George  have  inspired  very  interestiijg^^isQtrlftoa  Gaussian  elimination  and 
graph  algorithms,  see  Lipton,  ^osefiandiTarjan  [21],  Gilbert  [15]  and 
Rose  and  Whitten  [30].  In  our  stiidy,  Wfi^vMse  the  language  of  block 
Gaussian  elimination  extensively,  but  we  focus  exclusively  on  iterative 
methods  which  do  not  require  that  the  last  stage  of  the  Gaussian 
factorization  process  be  carried  out.  In  each  iteration  it  is  required 
that  all  the  subproblems  be  solved  with  somewhat  special  data.  Much  of 
this  work  can  be  carried  out  asynchronously  and  in  parallel. 

In  the  second  section,  we  consider  the  basic  structure  of  the 
linear  systems  of  equations  arising  in  a  finite  element  discretization. 
To  simplify  our  notations,  we  restrict  our  discussion,  throughout  this 
paper,  to  the  case  when  a  region  is  partitioned  Into  two  substructures. 
The  extension  of  our  results  to  cases  with  a  number  of  nonintersecting 
cuts  is  immediate. 
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In  the  third  section,  we  discuss  preconditioned  conjugate  gradient 
methods  in  general  and  introduce  three  iterative  methods.  We  show 
informally  that  two  of  these  can  be  expected  to  be  optimal,  i.e.  that 
the  number  of  iteratiottslf,;  r«ijulTesd'j  far  a  given  tolerance  will  remain 
constant  when  the  f_iniJ;ec«leriient rriiDdel^ls  refined. 

In  the  fourth  section,  a  . dedall^d  c.nalysis,  including  specific 
bounds  on  the  spectrumi of ithei^ltecatian  operators,  is  given  for  a  model 
problem,  namely  Poiaspn's  equaCionqen  L-  <ind  T-shaped  regions.  In  the 
fifth  section,  this  analysis>'  ts-  extended,  using  quite  different 
mathematical  tools,,  to  general rplane  regions  and  arbitrary  conforming 
Lagrangian  finite  elecent  cijaippDOxiifiationi;  of  second  order  linear 
elliptic  problems.    ,  -  ji^  Lirs^dorc  "  t:  a> .: :.  •:-  ;. 

Earlier  work  on  algorithms  of  this  general  nature  is  reported  in 
Concus,  Golub  and  O'Leary  [7],  who  gave  an  algorithm  and  numerical 
results  for  Poisson's  equation  on  .T^-'^stoffed  regions.  Two  additional 
preconditioners  were  recently  idirtroduced  in  Golub  and  Meyers  [16]. 
Methods  inspired  by  Schwarz's  alternating  principle  and  by  control 
theory  are  considered  in  HihUj  Glowinski,  and  P^riaux  [9].  Some  of 
these  methods  are  described  and  analyzed  in  section  6.  We  also  note 
that  we  have  learned  much  from  the  papers  by  Dryja  [11,12]  and  from  his 
October  1981  visit  with  the  second  author.  At  that  time,  we  were  first 
introduced  to  the  method  which  has  proven  most  successful  in  our 
numerical  experiments. 

The  authors  are  currently  actively  extending  their  work,  to  more 
general  problems.  We  have  begun  a  series  of  experiments  with  actual 
engineering  structures  and  have  also  made  progress  towards  the 
extension  of   the   theory   to   a  larger  class  of  elliptic  problems,  to 
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cases  in  dimensions  higher  than  two  and  to  non-Lagrangian  finite 
elements.  The  use  of  preconditioners  rather  than  exact  solvers  for  the 
subproblems  has  recently  been  considered  by  Bramble,  Pasciak.  and  Schatz 
[5].  Their  work  has  been  re  in  t-erprei-eii  In- blotk  Gaussian  elimination 
terms  in  Widlund  [34],  where  a  reiBtiobs.hip  --between^  their  algorithms 
and  those  of  this  paper  is  also  estsfellsbed.  f?  ,nox.-o9c  < 

We  note  that  both  from  an  algoEithtaiC-aad  analytic  point  of  view, 
our  methods  have  much  in  common  with  capacitance  -matrix  methods,  see 
Bjjirstad  [2,3],  Dryja  [10],  O'Leary  aad  Wldland.  [24j25] ,  Proskurowski 
and  Widlund  [26,27]  and  the  references  gluten- tnithose  papers.  See  also 
Widlund  [33],  in  which  related  :^£*tg<arithms  :f  or  mixed  finite  element 
methods  for  the  biharmonic  and  Stokes'  problems  are  also  discussed. 
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2.    SUBSTRUCTURED  FINITE  ELEMENT   PROBLEMS   AND  BLOCK  FORMS   OF   THE 
STIFFNESS  MATRICES 

To  simplify  the  discussion,  we  confine  ourselves  to  problems 
defined  on  a  connected  region,  ^,  which  is  a  union  of  J^  ^  ,  ^2  ^^'^  ^3 
and  to  the  Dirichlet  case.  Here  Ji  1  and  ^2  ^^^  plane,  bounded, 
nonintersecting  regions  and  To  the  intersection  of  their  closures.  The 
boundaries  of  fi  j^  and  ^2  ^^^  ^  1  ^  '^  3  ^^^  ^2  ^  ^3  »  respectively,  and  the 
boundary  of  ^   is  T  j^  u  F  2  .   We  assume   that   these   subregions   are 
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curvilinear  polygons,  i.e.  in  particular  they  are  Lipschitz  regions; 
see  section  5  or  Grisvard  [17].  A  linear,  second  order,  positive 
definite,  self adjoint  elliptic  operator  is  defined  on  ft .  Its  symmetric 
bilinear  form  is  denoted  by  aQ(u,v).  A  simple  example  is  given  by  the 
Laplace  operator  with  a  homogeiieo|as-.  Dixtchlet  condition  for  which 

.  -ajjCujv)  =?  />^o^^v  dx  ,   u,v  e  HqC^  )  . 

n 

Here  hA(JJ  )  is  the  subspace  !<afe-*el6m[e"nts  with  zero  boundary  values  of  the 
Sobolev  space  H^(n)  of  .^squaire":  integrable  functions  with  square 
integrable  first  der;i.v^tiv#3i.  cJ  '^riangulations  of  ftj  and  ft  2  ^^^ 
introduced  in  such  a  way  that  tjhe- node^  oil  F  ^  coincide  and  T^  follows 
element  boundaries.  We  assume  that  each  degree  of  freedom  of  the 
finite  element  subspaces  is  asso<iiateed.'  with  3  node  and  a  basis  function 
in  the  finite  element  space  and  that  the  support  of  any  basis  function 
coincides  with  the  triangles  to  which  ±x.s  node  belongs.  An  element  of 
the  stiffness  matrix  has  the  form  aQ  (t  i,i}> -j),  where  4)^  and  <^  i  are  basis 
functions,  and  it  therefore  vanishes  unless  the  two  nodes  belong  to  a 
common  triangle. 

In  the  case  of  a  general  curved  boundary,  the  region  ft  will 
normally  be  approximated  by  a  union  of  straight  or  isoparametric 
elements.  We  will  discuss  neither  the  modifications  of  the  bilinear 
form  which  is  requied  in  such  a  case  nor  the  related  issues  of 
numerical  quadrature  and  nonconforming  finite  elements;  see  e.g. 
Ciarlet  [6].  For  our  purposes  uniform  Vj^-ellipticity  of  Che  bilinear 
form  or  the  modified  bilinear   form  will   be   sufficient.  We  will 
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simllarly  assume  that   the   triangulation  of   the   region  Is   regular 
enough. 

It  is  easy  to  see  from  the  definition  of  the  bilinear  form  that 


J-c  ^.    .(V,u) 


ajj  (u ,  V  )  =  aj^  ^'(U^  V^ "  +^  ^ M.ui  V^^ 


iV.in^.''    r.   ;■>:.- 


(2.1) 


and  that  therefore  the  stiffnessL  macrfx  .  ofv  a  problem  on  ^  can  be 
constructed  from  those  of  J^  j^  and  ^2-  '^^  same  relation  holds  for  any 
pair  of  nonoverlapping  subregions. r 3  ^Eb±s ^f act  is  frequently  used  in 
practice  to  construct  stiffness  matrices  frdra'the  stiffness  matrices  of 
smaller  substructures.  A  special. case  is  -tfhe  common  process  in  which 
the  stiffness  matrix  is  assembled  f«daii''  t^fe'  contributions  from  the 
individual  elements.        :.  '.-dl3      ::b:'-i      an  r.--.. 

The  stiffness  matrix  of  the  entice  ln:-i>bl em  is  of  the  form, 

rr3  3i=rr  h-  . 


K  = 


K 


11 


o:  .-^--K- 


13 


K22-'  k:23 


^ 


K 


13   ^23 


K33  j 


where  Kj^j^  ,  represents  couplings  between  the  pairs  of  nodes  in  ^  j^  ,  Kj^^ 
couplings  between  the  pairs  belonging  to  H .  and  T  o  respectively  and  Koo 
couplings  between  nodes  on  I"  3  ,  etc.  By  (2.1)  this  matrix  can  be 
assembled  from  the  stiffness  matrices  corresponding  to  ii  1  and  ^ 2 
respectively, 


:(i)  = 


K^^   1C^3 


K 


13 


4P 


and 


■(2)  = 


^22  '^23 


I.  K.23 


4i) 


(2.2) 
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where  ^23       ^"^  ^^3     contain  only  the  contributions  from  the  Integrals 
over  fJ  j^  and  ^2  respectively.   By  (2.1),  we  see  that 


K33  =  41)  +  42)  . 


(2.3) 


We  note  that  the  degraes  of  f  re'6dom<  hSve;  been  partitioned  into  three 
sets  of  which  the  third,  the  separator  set,  corresponds  to  the  nodes  of 
To  .  From  the  point  of  view  of  graph  theory,  the  undirected  graph  of  K 
becomes  disconnected  into  two  components  if  the  nodes  of  the  separator 
set  and  their  incident  "edges "  dire  ^-"^femovedo  Since  conforming  finite 
elements  are  used,  i.e«  the  finite  element  space  V,  C  H  ,  it  also 
follows  from  the  assumptions  on  .aQ(u,v)  that  K  is  positive  definite, 
symmetric  and  as  a  consequence  so  are  Kj^j  ,  K22  and  K-j^  . 

The  matrix  K^^^,  defined  in  (2.2),  is  the  stiffness  matrix  of  the 
elliptic  problem  on  Q ,  with  a  natural  boundary  condition  on  T^  and  a 
Dirichlet  condition  on  F ,  ^  i.e.   the  problem 


^,(%'^h)  =  /  f  v^,  dx  +/  g^v^  ds  ,  V  v^,  e  V^  n  Hj(J^i,ri), 
J-  o  r 


(2.4) 


u^  £  V^  ,  Y  lUh  =  go  . 


Here  Y  j^u^  is  the  trace  on  F  j^  ,  i.e.  the  restriction  of  Uj^  to  F  j^  ,  and 
H^(f],,r,)  the  subspace  of  H  (f2]^)  with  vanishing  trace  on  F  j^  .  In  this 
case  Xo  is  a  vector  of  unknowns  and  the  linear  system  is 


<"  K^^   i<:^3 


M3  ^ 


wl 


^1 
^3 


I  b3 


(2.5) 
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We  note  that  the  vector  bj^  vanishes  if  f  and  gjj  are  zero  in  which  case 
the  vector  bn  represents  the  Neumann  data  on  F  2  . 

Written  in  variational   form,   a  Dirichlet  problem  on  fJ  j^  has  the 


form 


(2.6) 


^h  ^  \   ,     Tu^^  =:  §Dr»,.ng  ir  ^^i  iv  ic 

....;: nei.ijCDJOC    "Vj    orni    ;■■. 
where  Yu^   is    the   trace  on  ^  i^   T  2.  _  ^.lUymatrix  f  arm,  =.this   problem  can  be 
written  as  .       ,.,r.     -.-•-'-•-      «-'i         ,  =.    j 


"^11        ^^13    '' 


0 


^1 

^3    J 


r 


I-':?     '.'JfV).;^ 


■  ^  .-Ji^  riL.    c .  oi. 


t;     •   fti/r 


(2.7) 


"1 
lb3    J 

We  consider  the  linear  system  of  algebraic  equations  of  the  form 

.   I'or'  ■.  'I  :--r 


Kx  = 


^11 

0 

1  ■-  •■:  r 

^1 

0 

K22 

•^23 

^2 

'^13 

^^23 

K33 

^3 

By  using  block-Gaussian  elimination,  we  can  reduce  this  system  to  the 
positive  definite,  symmetric  system. 


Sx^  =  (^33~  ''^13  1^11  ^13  ~  ^23  ^^22  '^23'^3 


^3  ~  ^^13  '^11^1  ~  ^^23  ^^22  ''2  ~  ''3 


(2.8) 
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It   is   common  practice   to  complete  the  process  by  solving  (2.8)  by  a 
direct  method. 

The  right  hand  side  bo  can  be  obtained  at  the  expense  of  solving 
the  two  subproblems  on  iij  atld'Or^--',  mifltiT)lying  the  resulting  vectors  by 
the  sparse  matrices  RJ^  >  and  K2'3  -  respectively  and  subtracting  the 
resulting  vectors  from  bo.  From  now  on,  we  will  consider  only  the  case 
when  b,  and  b2  are  seTO."'  Such  a  red'uction  can  of  course  easily  be 
accomplished. 

The  matrix  S,  which  is  a"%b-ciilled  Schur  complement,  see  Cottle 
[8],  can  be  expensive  to  compute  and  store.  However,  we  notice  that  Sy 
can  be  computed  for  a  given  vector  y  at  the  expense  of  solving  the  two 
subproblems  with  the  sparse  right  hand  sides  Kj^^y  ar^d  ^^23^  » 
respectively,  and  certain  sparse  featirix  and  vector  operations.  In  the 
next  section,  we  will  develop  iter&ti^^e  methods  which  only  require  S  in 
terms  of  such  matrix-vector  products.'  <  - 

The  cost  of  computing  Sy"depends  primarily  on  the  efficiency  of 
the  solvers  for  the  subproblems.  It  should  also  be  noted  that  if  a 
Gaussian  elimination  method  is  used,  advantage  can  be  taken  of  the 
sparsity  of  the  vectors  Kj^^y  and  K23y.  Thus  when  the  lower  triangular 
systems  of  equations  are  solved,  the  computation  can  begin  with  the 
first  equation  which  has  a  nonzero  right  hand  side.  Similarly,  the 
solution  of  the  upper  triangular  systems  can  be  stopped  as  soon  as  all 
the  components  of  '^lI'^My  ^"*^  '^22'^23y'  necessary  for  computing 
^is^^ll^lSy)  and  K23(K~22'^93y ) »  ^^^^  ^^^"  found.  This  can  effectively 
reduce  the  size  of  the  triangular  systems  necessary  to  carry  out  the 
iteration  steps.  It  is  thus  particularly  advantageous  if  all  the 
variables  at  nodes  adjacent  to  F ^  are  ordered   late.    It   should  be 
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noted,  however,  that  such  a  constraint  may  be  hard  to  impose  on 
existing  software  or  that  it  may  lead  to  an  increase  in  the  time  and 
space  required  to  factor  K^^  and  K22  into  their  triangular  factors. 

We  will  also  need  the  5chur  cong)l,ement,s  .  with  respect  to  the 
matrices  K^^)  and  K^'^)  defined  \i^^,JJ^2\.:y'Ihey^  are,.. 


r!V  ,r!-:  won  •:;'ji  5 


■  a  cr-: 


s(l)  =  K^l)  -  k'[3K];[Ki3  and^.  s/2.)  ^   K^p  -  KI3K22K23,     (2.9) 


Using  (2.3),  (2.8)  and  (2.9),  we  f ipd. t)ia.t  ^i  r'  1  v 


::?.    .  : :     ■-   jjqr 


S  =  S(l),^3(^>:..v  :: 


(2.10) 


J.'  er 


:(1) 


The  mappings  S  and  S'^^''  play  an  imp9?'pant.  role  in  what  follows.  The 
vector  S'^-'y  can  be  computed  by  ,iSolvii»g. the  Dirichlet  problem  (2.7)  and 
then  applying  the  matrix  of  (2.5),  which  corresponds  to  a  Neumann  case, 
to  the  solution  vector.  This  can.^be  seen,  by  a  straightforward 
computation: 


K^3  K3^   _ 


•^11  '^n 


0 


-ir  0  ^ 
y 


°   1 

[s(l)yj 


(2.11) 


This  map  thus  takes  the  discrete  Dirichlet  data  y  on  T  3  into  the 
discrete  Neumann  data  s(^)y  on  the  same  set.  For  a  discussion  of  the 
continuous  case,  see  the  next  section. 

In  section  4,  we  return  to  a  discussion  of  the  structure  of  this 
problem  in  the  special  case  when  Q,  and  ^2  ^^^  rectangles  and  P 3  is  an 
interval.   We  have  carried  out   numerical   experiments   for   this 
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partlcular  geometry  using  uniform  meshes   of   right   triangles   and 
piecewise  linear  continuous  finite  element  functions. 

3.    CONJUGATE  GRADIENT  ALGORITHMS  FOR  SUBSTRUCTURED  PROBLEMS  AND  AN 
INFORMAL  THEORY..  .[-M.sr  ii  sr'.-'  .'I   i:  f  ' 

In  this  section,  we  will  introduce  three  iterative  methods  of 
conjugate  gradient  type.  The^^neral  theory  of  such  methods  is  quite 
well  known  and  it  will  therefore  be  discussed  only  briefly;  see  e.g. 
Concus,  Golub  and  O'Leary  [7],  Hestenes  [19]  or  Luenberger  [22]. 

Let  Ax  =  b  be  a  linear  system  of  algebraic  equations  with  a 
positive  definite,  symmetric:  matrix  A.  Let  x'  ^^  be  an  initial  guess  and 
r(0)  =  b  -  Ax^'^)  be  the  initial  residual.  The  k-th  iterate  in  the 
standard  conjugate  gradient  method,  x^*^',  can  then  be  characterized  as 
the  minimizing  element  for  the  variational  problem 

min(l/2)ytAy  -  y'^^b  , 
j 
where  y-x^^^  varies  in  the  linear  space  spanned  by   r^  ',   Ar^  -',   ..., 
^k-lj.(0)^    By  expanding  in  eigenvectors  of  A,  it  can  be  established 
that 

(x(k)_x)T  A(x(^)-x)/(x(0)-x)T  A(x(0)-x) 
is  bounded  from  above  by, 


min     max   (l-Xp(A))^  ,  (3.1) 

pGP^_l  Xea(A) 
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see  Luenberger  [22],  Here  x  is  the  exact  solution,  P^_]^  the  space  of 
all  polynomials  of  degree  k-1  and  a (A)  the  spectrum  of  A.  This  bound 
can  be  used  to  establish  that  the  convergence  is  rapid  if  A  is  well 
conditioned  and  that  the  rate- of  convergence  can  be  bounded  uniformly 
for  entire  families  of  operators  if  all  the  eigenvalues  fall  in  a  fixed 
interval.  .-..-lani   Ji.^v  ''^■'       ii-:ol2i,t-- 

Preconditioned  conjugate  .'gratiiient  methods  have  been  studied 
extensively  in  recent  years.  The  i-djea^^ses  back  to  the  mid-fifties; 
see  Hestenes  [19].  Let  Aq  bed aaother  positive  definite,  symmetric 
operator  for  which  it  is  feasible.ntpc solve:  iaaxiliary  systems  of  the 
form  Aqy  =  c  repeatedly  for  differfentrioright:  hand' sides.  In  one  of  the 
versions  of  the  method,  the  original- problem. Ax  =  b  is  transformed  into 
AA^^y  =  b.  The  iterate  y^"^' 'is,gought  as.  the  sum  of  an  initial  guess 
y(°)  and  a  linear  combination  pf.F^°],  Mg^r^^^  ...,  (AAgb'^'^r^.  If 
an  appropriate  inner  product  is  used,  a  convenient  recursion  formula 
results,  see  e.g.  Proskurowski  and-Widlund  [27],  Each  step  of  this 
algorithm  requires  the  solution  of  an  auxiliary  linear  system.  It  is 
important  to  note  that  the  estimate  (3.1)  still  holds,  but  that  now  the 
eigenvalues  of  AAq\  i.e.  those  of  the  symmetric  generalized  eigenalue 
problem  Ac})  =  XAq(j)  ,  are  of  relevance  rather  than  those  of  A.  It  is  also 
worth  noting  that  the  estimate  (3.1)  can  be  used  to  show  paticularly 
rapid  convergence  if  the  eigenvalues  are  clustered. 

In  problems  such  as  those  considered  in  this  paper,  it  is 
convenient  to  use  a  version  of  the  algorithm  in  which  the  operator  A 
only  appears  in  a  subroutine  which  provides  the  product  of  A-Ag  times  a 
vector. 

For   the  problem  at  hand,   we   first   consider   the  solution  of 
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equation  (2.8)  without  preconditioning.  From  (2.10),  we  see  that  Sy 
S^^^y  +  S^^\ .  It  is  therefore  at  least  plausible  that  S  will  be  ill 
conditioned  if  S^^^  and  S^^^  are.  As  shown  in  section  2,  S^'^^ 
represents  a  Dirichlet -Neumann  map  and  therefore  involves  a  loss  of  a 
derivative  in  I-^C^i)'  Such' 'a  "map  will  have  a  spectral  condition  number 
proportional  to  -the  number  of  nodes  on  F n.  In  order  to  clarify  this 
point,  we  consider  the  continuous  case,  leaving  the  details  on  the 
finite  element  case  to  sections-  4  and  5. 

Thus  consider  two  harmonic  functions  u,  and  U2  defined  on  fi  ^  and 
^2  respectively.  These  functiohs  vanish  on  F  j^  and  F2  and  have  the  same 
trace  on  F ^  .  They  can"  thetef&tfe  be  combined  to  form  u  e  HQ(n  ).  This 
function  satisfies     -  '^'     ■;  Jjri '--='■  ^' - 

;  -r  i  "ft.i.c  sill-  '  -■ 

an(u,v)  =/  Vu«^v'dx'='f(v)  ,   V  ve  hJ(Q  )  ,      (3.2) 

Q 

1-:.-'  r;n-.?' 

where  the  linear  functional  f  has  its  support  on  F  ^  .   It   is  easy   to 
show  that 


a^(u,v)  =  /  Vu^'Vv  dx  +  /  Vu2*Vv  dx  = 
"1  "2 

r  9^1    ,     r  5^2    ^     r  .3u,   . 

F3        r3        1  3 

where  v   is   the  normal  outward  with  respect  to  fJ  j^  .   We  can  therefore 
rewrite  equation  (3.2)  as 


(-9  u, 


[,^—]    =  f  ,   on  F  3  , 


-15- 

9u  ^^1        ^"2 

where  [^ — ]  corresponds  to  Sy.   Similarly  -r —  and  -  -r —  correspond  to 

S^^^y  and   S^'^^y   respectively.   Following  Lions  and  Magenes  [20]  and 

Grisvard  [20],  we   see   that   for  ue  H^Cil),  y -^u.     £   H^^^(r3),   where 

1 /2  1/2 

^00^^3^  is  a  subspace  of  H  '  (T^);    spe  .-section  5. 


It   then  follows,  from  a  standard  .variational  argximent ,  that 


9  U]^ 
IT 


^2      3  u  1/2 

_ —  and  [^— ]  belong  to  the  dual  space, , -pf ^ flQ^'rCr -j)  ;  i.e.   a  derivative 

is  lost  in  comparison  with  Y-j".    -rr:  stoi.x  l::nc ;  -"i  :  :=.  . 

In  view  of  what  we  have  just  learaed-^r  it- Is^  natural  to  try  to  find 
a  preconditioner  which  also  involves,  tjie  r  loss  of  a  derivative  in 
L2(ro).  A  natural  choice  would  be  a  tangential  derivative  but  that  is 
not  a  symmetric  operator.  Instead  we  can  use  the  square  root  of  the 
negative  of  a  discretization  of  the  Laplacian  on  F  ^  ,  Such  a  method  is 
practical  at  least  for  problems  in  the  plane  and  has  been  tested;  see 
sections  4  and  7  for  details.   We  denote  this  operator  by  J. 

In  our  experience,  an  even  better  method  involves  the  solution  of 
a  system, 


c  t ' 
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This  solution  can  be  obtained  by  solving  equation  (2.5)  with  the   right 
hand   side   of   (0,y)    and  then  the  discrete  Dirichlet  problem  on  J^  2  . 


using  x-j  as  data,  cf.   (2.7).   The   relevant  mapping   is  now  SS 
since 


(1) 


-1 
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■  Hi 

0 

K.^3 

"^11 

0 

•^22 

K23 

0 

K 

'^23 

K33     j 
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.22  1^23 


4\'i 


^  -1 

'  0' 

0 

= 

.  y  . 

ss 


0 
0 

(1) 


-1 


For  an  analysis  of  the  spectrum  of  SS^^-*   ,  see  sections  4  and  5. 

We  note  that  we  need  not  construct  any  auxiliary  operator  when 
this  preconditioner  is  used.  It  is  also  interesting  to  note  that  if  a 
symmetric  region  is  cut  in  half,  and  treated  fully  symmetrically,  then 
S  =  2S(^)  and  the  conjugate  gradient  iteration  converges  in  one  step. 

4.    ANALYSIS  OF  A  MODEL  PROBLEM. 

In  this  section,  we  consider  a  model  problem  where  fi,  and  ^2  ^^^ 
rectangles  and  ^  is  L-  or  T-shaped.  The  domain  is  triangulated  making 
right  triangles  of  equal  size.  The  number  of  interior  nodes  in  the 
horizontal  and  vertical  directions  are  q  and  r  respectively  for  ^  ^  and 
m  and  n  for  ^^2*  There  are  p  internal  nodes  in  the  horizontal  direction 
before  a  node  in  f2  2  aligns  with  the  first  vertical  column  of  nodes  in 
^1  .  We  assume  that  p  +  q  £  m.  With  piecewise,  linear  continuous 
elements  and  Poisson's  equation,  we  obtain 

•^11  =  (Ir  ®  Rq)  +  (^r  «  ^q  >  ' 
^^22  =  (^n  ®  \,)  +  (^n  ®  ^m)  , 


where  R  =  tridiag(-l , 2, -1) .  Furthermore 


^13  =  (O.-Iq)'  ' 
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T 


K23  =  (0,-Iq,0) 


and 


^3^  =4P  =^  t.,l^V:--..o:  .r:  I 


We  will  also  use  the  notation  >J_,  =  Rq  -•   This  is  a  positive  definite 
symmetric  matrix. 

The  orthogonal  matrix  Q„  defined  by 


(Q„)i.    =  (2/(q  +  l))^/2   sinCljTr/(q  +  l)] 


^q'ij 


provides   the  eigenvectors   of^  R  and  J^  .and  QgRnQa  ~  ■'^a  ~  diag(Xi^^) 
with 


..ni*r  .'■.•'i 


X^q)  =  4  sin2(jiT/2(q  +  l)],  j  =  1,2, ...,q. 


The  eigenvectors  of  K,.  and  K22  can  be  constructed  as  tensor  products. 
Thus 

(Q^  8  Qq)  K^i(Q^  8  Qq)  =  (I^  8  A^)  +  (A  ^  «  I^)  =  ^-qj- 

I       ' 
and  similarly  ^ 

(Qn  ®  Qn,)  K22(Qn  «  Qm)  =  ^^n  «  ^m)  +  ^^n  «  In,)  =  f'mn  • 
With  the  Schur  complement  S^^)  defined  as  in  (2.9),  we  find  that 
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QqS^^\    =   Iq    +    1/2  Aq    -   Q^k[^(Q^  «    Q^ )   il^'^CQ^  «    Q^^i^Q^ 


r 
=   Iq    +    1/2  Aq    -    (2/(r+l))     Z      sin^f  tar /(r+1 ))     (Aq+  X^"^)!)"! 

k=l 


=  I    +  1/2  A    -  D      .■  y-^  ^ 


Lemma  4.1.   The  elements   off  ttie  diagonal  matrix  D   ,  defined  above, 
are  given  by 

d(q.O  =a(0[l  -  (a(q))2r)/(l  -  (a(q))2(r+l  )^  ^   ^  =  1.2.... .q. 


where 


i(q)  =  1  +'x(<^)/2  -  (X(l)(l  +  x(<l)/4))  1/2.    (4.2) 


Proof:    The   matrix  D        is   diagonal   and   its   elements   are   given  by. 

r 
(l(q.r)   =   2/(r+l)  Z      sin2(k:r/(r+l)) /(  4  sin^C  jir /2(q+l  )  )   +  4sin2(lcir /2(r+l  )  )) 
k=l 

r 
=   (2a(q)/(r+l))  Z      sin2(kTr/(r+l))/(l-2a('1>cos(lciT/(r+l)   +   (aij^h^). 
-■  k=l  ^ 


This  sum  can  be  evaluated  by  using  Poisson's  summation  formula.   Thus 


d(<l'r)  =  (2aj^'l)/(r+l))[(l/2)f(0)+f(TT/(r+l))+f(2^/(r+l))+...  +  (l/2)f(7r  ): 


(2a(q)/Tr)  [Fq  +  2Z   F2i,(,+  i)]. 
k=l 


where 
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2 
f(x)  =  sin^x  /  (1  -  Za^'l^  cos  x  +  a^^^  ) 


and 


IT 

F  =  f  f(x)  cos  mx  dx  .    r 
0 


By  using  Poisson's  integration  formula,  w,er  otvtain 


Fq  =W2 

F^  =  Tra(^V4 


F^  =  -  Tra(^)   ~   (1  -  a(^>  )  /  4  ,   m  =  2,3,, 


Therefore 


d(q'^)  =  a(q)(l  -  (a(q)  ^  -  1)  E   (a(^))2('^+l)k) 
J        J         J         k=l   ^ 


=  a(0(l  -  (a('^))2r)/(l  -  (a^^h^^'^'h. 


The  parameter  ai'l',  defined  by  (4.2),  decreases  monotonically  from 
1  to  3  -  2/2"  as  X\'^^  varies  from  0  to  4.  The  ratio  d\^'^^/a^^^  is  also 
monotone  as  a  function  of  a>^'  and 


(r/(r+l))  a^^q)  <  dj^'^^  <  aj^^  <  1  .  (4.3) 


Denote   by  E     the   Ttf<q   matrix 

m 


E^  =    (0,    Iq    ,    0)T 
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where  the  leading  zero  block  is  p^q.      Similarly  as  in  our   computation 
of  the  eigenvalues  of  S^  •',  we  obtain, 


\^^^\    =  Tq  -^  1/2  Aq  "  Qq43(Qn  ®  Qm^  Vi^Qn  ®  Qm)  ^ZS^q 


Iq  +  1/2  Aq  -  Q^K%^mn%^m% 


T  '  -  - .   '' 


r  f 


Lemma  4.2.   For  all  x 


^^^m^mEm^  <  ^^-^q^  « 


Proof:  We  will  show  that  a  principal  minor  of  the  square  root  of  a 
positive  definite  matrix  is  dominated  by  the  square  root  of  the 
corresponding  minor,  from  which  Lemma  4.2  follows. 

Let 


B  =  A 


1/2  = 


^11   1^12  ^ 
^12  ^22  J 


Then   A 


11 


J^p2 


^11   "*"   ^12^12   ^^^        •'■'■   ^^        easy   to   see   that 
J-^1/2, 


x^B^^x  <  x^(Bf^+B^2Bl2)   ^• 

We  will  now  give  upper  and  lower  bounds  for  the  quadratic  forms 
x^S^l^x  and  x^S^2)x,  which  can  be  used  to  obtain  close  bounds  on  the 
performance  of  the  conjugate  gradient  methods.  We  first  note  that, 
with  y  =  QgX,  XX  =  1, 


x^S^Dx  =  y^QqS^DQqy 
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=  y^dq  +  1/2  Aq  -  Dq^)y  . 


By  a  simple  calculation,  we  see  that, 

An  upper  bound  is  similarly  obtained  by  using  (4.3): 


Js^^x   <  yT(i^+   1/2  A^)y-   (r/(r+l))   y^(Iq+    1/2  A^   -  ^lA\^^)y 

<  (l/(r+l))y'^(Iq+    1/2  ^^)y+  /I  y^Aj/^y  (4.5) 

<  (3/(r+l))   +  /I  x^JqX   . 


We  will   also   use   the   fact   that 

•  ■: .,    1-  ■  .■  r>r  ^  ■?  ? 

•  .  '  i. '  .:  r.  .  ii 
x'^s(2)x   >   0 


and  that  by  (4.3)  and  Lemma  4.2, 


xTs(2)x  <  yT(I   +  1/2  A   -(n/(n+l ))(!,+  1/2  A_  -/2  QqE^QnA^  ^Qn^E^Qn )  )y 


<  3/(n+l)  +  /2  x'^E^J  E 


m  m~m^ 


<_  3/(n+l)  +  /2  x'^J  X 


To  obtain  an  estimate  for  the  preconditioner  J   ,  we  note  that 


1  =  x^Jqx/x'^JqX  <_  x'^s(l)x/xTj^x  <  x''^(s(^)+s(2))x/xTj^x 
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<_  (3/(r+l)  +  3/(n+l)  +  Z/I  x'^JqX)/x'^JqX 

<  (3(q+l)/Tr)  (l/(r+l)  +  l/(n+l))  +  2/T. 

If  S^^^  is  used  as  a  preconditioner ,  we  obtain  the  estimate 

■     '■''■  ^ 
1  <  x'^^CS^l)  +  S(2))x  /  x^sC^x 


<_  1  +  (3/(n+l)   +  /2  x'^JqX)/x^JqX 


<  1  +  3(q+l)/Tr  (n+1)  +  /2  . 

The  upper  bounds  go  to  infinity  with  (q+l)/(n+l).  Numerical 
experiments  indicate  that  a  more  careful  analysis  might  eliminate  the 
terms  containing  that  factor.  Similarly  experiments  indicate  that  the 
true  upper  bound  for  x^(s(^)+  s(^^)x/x^s( ^ ^x  is  2.  In  all  cases  which 
we  have  tried,  the  lower  bounds  have  exceeded  1.6.  See  further 
section  7. 

5.    BOUNDS  FOR  THE  SPECTRA  IN  THE  GENERAL  CASE. 

In  this  section,  we  will  formulate  and  rigorously  prove  our 
results  on  the  spectra  of  the  operators  S,  SJ~^  and  S(S^  ')"  which 
correspond  to  three  methods  introduced  in  section  3.  Almost  the  entire 
section  will  be  devoted  to  the  proof  of 

Theorem  5. I.  The  condition  number  of  S  grows  linearly  with  the 
number  of  degrees  of  freedom  associated  with  F  t.  Those  of  SJ  and 
S(S^^^)~^  are  uniformly  bounded  and  the  corresponding  methods  are 
therefore  optimal. 
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Our  results  will  be  established  for  conforming,  Lagrangian  finite 
element  approximations  of  Dirichlet  problems  for  self-adjoint,   second 
order  elliptic  problems  in  plane  regions, 

Lu   =  -     I a^j(x)  gJi-  +  aQ(x)u  =  f(x),   x  e  fi  , 

i,j   '^i         ^j 

u(x)  =  g(x),    X  e  r^  U  f  2  . 

The  operator  L  has   real,  sufficiently  smooth  coefficients,  a..(x)  = 
a^.(x)  and  the  bilinear  form        _ 

a(u,v)  =  /  I      a^j(x)  g_iL  g-L  +  ao(x)uv  dx 


satisfies 


f  IS   : ".  r 


j_  „^||2     <  a(u,u)  ,    Vu  e  hJc^  )  (5.1) 

and 

la(u,v)l   <    Cllull     ,  llvll     ,  .  (5.2) 


We  adopt  the  common  practice  of  using  C  to  denote  a  generic 
constant,  strictly  positive  and  independent  of  the  dimension  of  the 
finite  dimensional  problems  considered,  with  the  value  of  C  not 
necessarily  the  same  in  different  instances. 

By  assumption,  the  finite  element  problems  on  the  two  subregions 
are  solved  exactly.  We  therefore  concentrate  our  study  on  the  solution 
of  equation  (2.8),  or  what  is  the  same,  the  discrete  problem  with  zero 
boundary  conditions  on  T^   U  T2   and  a   right-hand   side   different   from 
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zero  on  r^  only.   Any  problem  can  be  reduced  to  such  a  problem  by 

solving  two  subproblems. 

Even  if  the  original  problem  has   a  very   smooth  boundary,   its 

partition  into  subregions  makes  it  necessary  to  consider  regions  with 

corners.   Our  main  technical   tools  will   be   borrowed   from  elliptic 

regularity  theory  which  is   very  well  developed   for   the   case  of 

sufficiently  smooth  boundaries;  see  e.g.   Lions  &  Magenes   [20].   Many 

regions  arising  in  practice  have  corners,  but  in  spite  of  this,  the 

extension  of  the  theory  to  such  cases  has   been   relatively  neglected; 

see  however  Grisvard  [17].   Followiiig  Grisvard  [17],  we  concentrate  our 

attention  on  curvilinear,  polygonal  regions.   The  boundary  T    of  such  a 

region  is    U    V  "^     where  F^  Is  the  closure  of  an  open  sufficiently 
finite   ^  -^  ' 

smooth  curve  FT  .  Locally  such  a  region  is  the  image  under  a  smooth 
mapping  of  a  region  with  a  corner  with  an  angle  equal  to  n  ,  it /2  or 
3ii/2.  Denote  by  C  the  class  of  I  times  continuously  dif f erentiable 
functions. 

Definition.   Let  ^  be  a  bounded  open  subset  of  R  .   Its  boundary  F 

Q 

is  a  curvilinear  polygon  of  class  C    ,  i    an  integer  >  1,  if  for  every  x 

2 
e  F   there   exists   a  neighborhood  V  of  x  in  R  and  a  mapping  ^    from  V 

o 
into  R^  such  that 

(i)   T|^(y)  =  (.^  Y^y),i>2(y))    is  injective, 

(ii)  t  and  t"^  e  C", 

(iii)  n  n  V  is  either  {y  t   Q    \   ^2(7)  <  0}  , 

{yea    I  t  i(y)  <  0  and  ^■yiy)    <  0}  or  { y  e  a  |  t  i(y)  <  0  or  1 2(7)  <  0^  • 

We  will  assume  that  I  is  large  enough  and  that  the  separator  curve 
Ft  is  sufficiently  smooth.  Our  arguments  easily  extend  to  the  case 
where  F^  is  only  piecewise  in  C  . 
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We  will  need  a  number  of   Sobolev  spaces   and   some  of   their 
properties.    By  L2(^),   we   denote   the   space  of   square  integrable 
functions  on  ^.      For  integer  m  >  0,  H™(Q )  is  the  subspace  of  L2(n  )   for 
which,  .   -  -  ;  . .  •   ,  ■   .. 

i  ' '::       e  Ac-:  ; 

"""»m,o.=(^    ^   |(^«u(-K)]2d^)l/2  <«   .         (5.3) 
"  ^^'^     n  |a|<m 

.  ■  .s   at'    ,  '■  •hi'-.i '  :■  •'.. 

Here  (i_)«  =  {±.)"H±-f\    |a  |  =  a  ^%- ti'^.   --^^' 

For   s=m+a    >0,    0   <  a    <1,   we   define  H^(Ji)   in  terms   of   the   norm, 

"""Hsr.^  =  ^"^"n-r.^"  ^  ^  ?  lU^"-^-')  -(^"-(y))  (x-y)-^^^>|2  dx  dy)  1/2. 

(5.4) 
iJiv  "  .:7O0 

We  will  also  need  the  seminorms  lul      ,  which  are  obtained  by 

'  'HS(n)' 

dropping  all  terms  with  derivatives  of  order  less  than  m  in  (5.3)  and 
the  first  term  in  (5.4).  We  can  of  course  also  define  Sobolev  spaces 
on  the  curves  T  ., 

For  Dirichlet  problems,  we  use  Ht(S^),  s  >  0,  which  is  a  subspace 
of  H^(n)  defined  as  the  closure  in  H^(i;^)  of  D(fi  ),  the  space  of  Cf° 
functions  with  compact  support  in  Q.  The  dual  space  of  H^(n  )  is 
denoted  by  H"^(Q  ). 

Let  Y^   be   Che   trace  operator  on  F.,   defined  for  continuous 

functions  as  the  restriction  of   the   function   to  T  ..    This   set   of 

operators   can  be  extended  to  an  operator  from  H^(fi)  onto  a  subspace  of 

n  H  '  (r.).  This  subspace  Is  characterized  in  detail  in  Theorem 
J       -J 
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1.5.2,3  in  Grisvard   [17].    Here  we  are  primarily  interested  in  the 
trace  of  functions  in  the  spaces  Hi(f2.,r.),  i  =   1,2,   defined  as   the 
subspaces   of  H-'-(Ji.)  with  zero  trace  on  T  ..   In  this  special  case  it 
follows  that  Y^  maps  E.q{Q^,T^)   onto 


=  (Uul|2       ^  llp~l/-ulu  (r,    0^^^<"}- 


Here  p  is  the  distance  to  the  boundary,  i.e.   to  the  end  points  of  F  .j. 

The  space  H^^'^Cr^)  is  strictly  contained  in  Hg'  ^3);  see  Lions 
and  Magenes  [20],  p.  66.  It  can  also  be  defined  by  interpolation. 
Thus,  .  .,,.  ..    .      ..^  ..  ,.  -.,-,,. 


Hj^2(r3)  =  [Hl(r3),  L2(r3)]i/2 


We  also  need 

H^/^(r3)  =  [Hj(r3),  L2(r3)]3/4  ; 

see  Lions  and  Magenes  [20],  pp.  64-66.  We  can  use  the  K-method,  see 
Lions  and  Magenes  [20],  pp.  98-99,  to  find  useful  formulas  for  the 
B.y^{V  2)  and  HQ/'^(r3)  norms  of  finite  element  functions.  Let  u  e 
H^(r3)  +L2(r3),  i.e.  u)  =a)Q+aj^,  where  o)  q  e  HqCF  3)  and  w  ^  e  L2(r  3). 
A  functional  K(t,u))  is  defined  by 


K(t.a))  =   inf   («t^o"^l,„    +  t2lla,il|2     )l/2. 


The  equivalent  norms  are  then  given  by  the  expression 
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1/2 
^^      .  +  I  t  ^"  "•"  ^(K(t.a)))^dtl 


^"•^"Lcr,)  +/  t-^^^^^\ut,i^))ht] 


with  6  =  1/2  and  3/4  respectively.  By  using  the  same  arguments  as  in 
Lions  and  Magenes  [20],  we  can  now  compute  the  ^"IAq  ^^3^  ^^'^  ^0  ^^3^ 
norms  of  o)^,  the  restriction  to  T^  of  a  finite  element  function  defined 
on  the  triangulation  of  Q.  .  -The  H^(r  3)  norm  of  ui^  equals 
(ap  ('^h>'*'h^^    '  ^''^s'^s  the  bilinear  form  ap  (u,v)  is  defined  by 


,    ,  du  dv 

ap    (u,v)   =  J    [-—  — 

3  p      ds   ds 

^  3 


P    (u,v)   =/    [JiJL+uv)   ds    »      u^v  e    Hj(r  3) 


For  a  given  finite  element  space,  we  can  "compute  the  stiffness  and  mass 
matrices  K^,  and  Mp  associated  with  this  Dirichlet  problem  on  T  3.  If 
the  coordinate  vector  of  o)^  with  respect  to  the  finite  element  basis  is 
a,  then  it  follows  from  the  computation  that 


"'^h"Hl/2,p    =a'^J«  . 

Hqo  (r3) 


where  J  =  (Mf  ^/^^Mf ^/2)l/2  ^^^   ^^^^ 


lla),ll  2        =  aT  j3/2 

^H3/^(r3) 


Let  Q  be  the  orthonormal  matrix  of  eigenvectors  of  J  and  A  be  the 
corresponding  diagonal  matrix  of  eigenvalues.  It  is  then  easy  to  see 
that 


aT  j2e  ^  =  3T  ^29 


where  a   =   Cp  .    We   also  note  that  the  pair  N^'-'^Q  and  A   solves  the 
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generalized  eigenvalue  problem  for  the  pair  of  matrices  ^    ,      ^  ,      i.e. 

In  what   follows,  we  also  need  to  know  that  the  trace  class  on  F o 
of  E^/^(a^)        Hj(fi^,r^),  i  =  1,2,  is  given  by  H^/^CF  3)  and  that 


r  -..^■■1 


'"'3-H3/V3/  ""'H5M(a^)  .   1-1.2.  (5.5, 


where  the  constant  depends  on  the  region  only. 

One  of  the  principal  results  in  elliptic  theory  is  sometimes 
called  the  shift  theorem.  For  Dirichlet's  problem  and  an  elliptic 
operator  L  of  order  2k  it  states  that  „.  ,,  • 


II  ull 


M  <    Cf  II  fll         ,  +  llunll       ,,  )      .  (5.6) 


Here  Lu  =  f  and  Uq  is  any  extension  of  the  Dirichlet   data  such  that 

u-Uq  e  Hq   (f^  ).   It  is  easy  to  establish  this  result  for  s  =  0  for  our 

family  of   second  order  elliptic   problems   by  using   a   standard 

variational  argument.   These   bounds   can   also  be  derived  for  a  wide 

range  of  values  of  s  for  problems  with  sufficiently  smooth  boundaries, 

but   this  is  no  longer  true  for  Lipschitz  regions.   However,  Necas  [23] 

has  shown  that  the  estimate  (5.6)  is   true   for   |s|   <   1/2.    In  the 

estimate,   the   term  II  u^H   ,,      can   be  replaced  by  a  suitable  trace 

0  HS+k(Q) 

term,  by  using  an  extension  theorem  which  states  that  there  exists  a 
continuous  operator  which  extends  any  element  in  the  proper  trace  class 
to  an  element  in  H^''''^(n);  see  Grisvard  [17]  or  Stein  [31]. 

Before  we  turn  to  the  proof  of  Theorem  5.1,  we  will  show  that  the 
continuous  operator  corresponding  to  S  S*^^'  is  bounded  in  an 
appropriate  norm.   For  simplicity,  we  will  use  the   same   notation   for 
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the  operators  in  the  continuous  case.   We  note  that  a  bound  for  s'^-^S" 
is  easy  to  obtain  since  S  =  S^^)  +  S^^^  and  S^^^  and  S^^^  are  positive 
definite  and  self adjoint. 

The  natural  domain  of  definition  fpr,(S^^b~^  is  (Hq^'^CI' 3) )',  the 
dual  of  the  trace  class  H^^-^CT^).   Consider  the  variational  problem 


"4 


^-^i  av 


a^.C'^l.v)  =  j   [    I    a^.  ^-—^  +   aou^v)  dx  =  /  g^v   ds  , 

1  ^^        J  dXj  dx.  .    j,^ 


1/2  •'  .  r  -    i       ■ 

Here  gjq  e  (Hqq  (r^))'  is  the  Neumann  data   given  on  T  2-        From  the 
ellipticity  assumption  (5.1),  we  obtain 

■■■  u  "  ■'-'      '     ■  C'  >''■  ^ 

r""i"Hira  ^  *"  ^,("i."i)  =/  Sn^i  ds  <  "gN"  1/2,.  ,,^  "^3"i"  i/2.r  , 

H  (S2^)      1  To  ^^00  ^^3)>        "00  *^^3) 


By  the  trace  theorem 


"3'"'hJ^V3)'"'"h'(«i)  •    »^=''i<«l''-l>' 


Therefore 


"^3"i"  i/2.r  ,'  ''"sn"  Hi/2rr  ^^'• 

Hqo  CI  3)         CHq5  {T  2)) 

By  the  characterization  of  the  trace  class  of  Hq(S^2>^2^>  ^^  ^^^  that 
Y3U]^,  extended  by  zero  on  F  £,  is  the  trace  of  some  U2  e  HQ(n2»I'2^»  ^^^^ 


-  h1(J^2)  «00  (^3) 
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By  using  the  shift  theorem  in  the  simple  case  of  s  =  0,   we  can  solve 
Lu2  =  0  on  5^  2  with  U2  -  U2  e  ^ki^  2^ '      ^^  note  that 

L(u2— U2)  =  -  Lu2  >    U2  ~  U2  G  Hg(fi2)  » 
and  that 

Thus  we  obtain  Uj^  e  HqCJ^j^jFj^)  and  U2  £  HQ(f2  2>^2^  which,  by 
construction,  have  the  same  trace  on  To.  Together  they  therefore  form 
u  e  H^(Q  )  with 

II  uii  ,   <  c  n  gxjii  1  /o  "^  ■ '  . 

The  operator  L  maps  HA(n  )  continuously  onto  H~  (J^  ).  In  this  case  Lu  is 
a  distribution  with  a  support  on  T^  and  can  therefore  be  regarded  as  an 
element  of  iE^(^'^ (T  ■^) ) ' .  The  boundedness  of  3(5^^^)"^  is  thus 
established. 

The  proof  in  the  discrete  case  proceeds  somewhat  differently, 
since  we  want  to  obtain  bounds  for  SJ~  as  well.  We  also  need  to 
establish  an  extension  theorem  in  the  finite  element  case.  The  proof 
of  Theorem  5.1  follows  from 


1  J  <  S*^^^  <  C  J  ,   i  =  1,2.  (5.7) 


To   establish   the   lower  bound,  we  can  argue  as  in  the  continuous  case 
and  consider 
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^3 


"i,h^  ^h  ^Hjcn^.r^)  . 


As  above  we  find 


"OO  <-^3->       1    '     '       P_^ 

By  using  the  definition  of   the  dual  nonn  and   the   formula  for 
"'^3"l,h"  l/2.r  ,'  "^  obtain 


'00  "^^  3 


2         *•  1-  *  T  T-l;! 


"^3-l,h"  l/2,p   <  C6^  J-6.  (5.8) 


where  the  components  of  6  are  J    gjq4)^ds,  where  ij)^  is  a  basis  function 

^3 
for  the  finite  element  space  on  F^.   The  coordinate  vector  for  Y3UJ  ^ 

is  (S^^O~  5  and  the  inequality  (5.8)  can  therefore  be  written  as 
5T(s(l))-lj(s(0)-l  6  <  C  5T  j-k 

from  which  the  lower  bound  of  (5.7)  follows  by  using  that  S^  ■^  and  J 
are  positive  definite  and  symmetric  matrices. 

The  main  difficulty  in  establishing  the  upper  bound  is  the  proof 
of  the  following  lemma.  ^ 

Lemma  5.1.  There  exists  an  extension  mapping  from  F  ^  to  !i^  1  which 
extends  a  finite  element  function  g,  on  Ft  to  an  element  vi^  e 
V^  n  H^(n  ^,r  j^)  such  that 


^'^'H'(n.,'^"^'>'Hj/^(r3)  •  ''■'' 
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Proof :   We  first  construct  v  e  HgC^j^.Tp  by  solving  the  continuous 
homogeneous  Dirichlet  problem  with   the   prescribed   boundary  values. 
Similarly,   we   define   v,   as   the   solution  of   the   finite  element 
approximation  to  the  same  problem.   By  the  triangle  inequality 


II  v.ll  1     <  ii  vii  1     +  II  vv,-vil  , 


The  solution  of  the  continuous  problem  can  be  estimated  as  follows, 


II  vll     ,  <    dig.  II     ,  /o 


and 

II  vll 


H5M(a^/    "'^h"H3/4^,^)    ; 


see  Grisvard  [17],  Chapter  1.5.2.   From  the  variational  formulation  of 
the  discrete  problem  and  (5.2),  we  know  that 


J-  2 

—  II  v,- vll     ,  <    ao    (vv,-v,Vi  -v) 

c       h      nh^i)       ^1     h     '  h     ^ 

/v,-vii   ,         iin.  v-vii    ,  , 


<■    C    II  V.  -vll 


where  H,  v  is  the  interpolant  of  v  in  V.  .  It  is  important  to  note  that 
H^v  -  V  =  0  on  r  j^  U  r  2.  We  use  the  following  lemma  to  complete  our 
proof  of  Lemma  5.1. 

Lemma  5.2.   The  interpolant  IT,  v  satisfies  the  bound 


in,  v-vll  ,     <  C  h^/^lvl  c/- 


-33- 

With  the  aid  of  Lemma  5.2,  we  conclude  that 


Hv.-vII  ,     <  C  h^/^|v|  c/A 


,1/^11  ..»  ...     <  c  hl/^(5Tj3/25)l/2, 


<  c  h^/^iig^n  ^1, 


It   is   easy   to   establish,   under  our  assumptions  on  the 
trlangulation,  that  |j|  <  Ch~^  and  therefore 

(5Tj3/25T)l/2  <  ^  ^^l^i^'^l^^^ 


=  C  h-1/^llgull, 


which  completes  the  proof  of  Lemma  5.1. 

Proof  of  Lemma  5.2:  The  proof  is  a  modification  of  the 
Bramble-Hi Ibert  lemma.  Let  F(u)  be  a  bounded  functional  on  H^'^(A), 
where  A  is  a  triangle,  such  that  F(p)  =  0  for  all  polynomials  of  degree 
one.   Then 


|F(u)|  <  C|u|  .,, 


Since  the  modification  of  the  standard  version  of  the  Bramble-Hilbert 
lemma  is  relatively  minor,  see  Ciarlet  [6],  details  will  not  be  given. 
We  only  note  that  the  result  is  obtained  by  scaling  and  adding  the 
contributions  from  all  the  triangles  in  ^,.  We  also  note  that  it  is 
easy  to  prove,  by  Sobolev's  lemma,  that  H,  is  a  bounded  operator  in 
H-''  since,  by  assumption,  we  work  with  Lagrangian  finite  elements  only 
and  thus  use  only  values  of  v  and  no  values  of  derivatives. 

With   the   help   of   Lemma  5.1,  we  can  reduce  a  discrete  Dirichlet 
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problem  to  one  with  homogeneous  boundary  values   for  which  a  simple 
variational  argument  is  sufficient.   We  can  thus  establish  the  estimate 
(5.10)  for  the  approximate  solution  u^  ^  of  Lu  =  0.  For  given  Dirichlet 
data  gQ  e  ^00^^  3^'  "^  ^^^   thus  find  uj^  j^  e  HgCf^j^.T^)  such  that, 

Translating  the  results  of  Section  3,  we  find  that 

^/^l,h'^)  =  /  S^^^gDV  ds  ,   V  V  e  H^(fJ^,rp, 
from  which  follows,  that 

'§d"  „i/2.r  )'  '  ^"Sd"  l/2.r  ^ 


^3 


lis(^)gjl   wo      <  CiigJI 


By  using  the   formulas   for  the  norms  of  finite  element  functions,  we 
obtain, 

If  ' 

s(i)  <  C  J  . 

The  proof  of  the  theorem  can  now  be  completed  since  the  bounds  for 

—  1         (\  "i~^ 
SJ  ^  and  SS^  -^   follow  immediately.   It  is  also  easily  seen  that  the 

condition  number  of  S  grows  like  that  of  J,   i.e.    linearly  with  the 

number  of  nodes  on  F  -,. 

6.   A  discussion  of  certain  other  methods. 


In  this  section,  we  will  describe  and  analyze  some  algorithms 
introduced  in  Concus,  Golub  and  O'Leary  [7],  Dihn,  Glowinski  and 
Periaux  [9]  and  Golub  and  Meyers  [16], 

In  a  well  known  paper  on  generalized  conjugate  gradient  methods, 
Concus,  Golub  and  O'Leary  [7]  considered,  among  other  things,  the  five 
point  difference  approximation  of  Poisson's  equation  on  a  T-shaped 
region,  i.e.  the  same  problem  as  in  sections  4  and  7  of  this  paper. 
To  describe  their  approach,  we  modify  our  notations  slightly.   We  split 
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off  the  q  mesh  points  in  JJ  2  >  immediately  adjacent  to  the  horizontal 
separator  set,  and  associate  the  index  4  with  this  new  set,.  The 
resulting  finite  difference  problem  is  then  of  the  form. 


•^11 

0 

K^3 

0 

^1 

■  bi  ^ 

0 

K22 

0 

K24 

X2 

f 

b2 

^h 

0 

K33 

K34 

""^ 

^3 

0 

■^24 

K.34 

5^44  . 

.      X4  J 

.  ^4  . 

(6.1) 


A  symmetric,  positive  definite  preconditioner  is  obtained  by  replacing 
the  K3^  and  K34  matrices  by  zero.  After  a  symmetric  permutation,  the 
preconditioner  is  seen  to  be  a  direct  sum  of  two  discrete  Dirichlet 
problems  on  the  rectangles  ^  ,  and  ^2  •  Without  loss  of  generality,  we 
can  set  the  subvectors  bi  and  b2  equal  to  zero,  since  we  can  reduce  the 
system  (6.1)  to  that  form  at  the  expense  of  solving  two  discrete 
problems  on  the  rectangular  subregions. 

The  rate  of  convergence  of   this   generalized  conjugate  gradient 
method  is  determined  by  the  eigenvalues  of  a  2q  x  2q  matrix. 


S  = 


-(I+(1/2)R  +s'^^))"^  1 


-(I  +  (1/2)R^+S<^^))"^ 


^34 


^^43 


Here  S^-'-''  and  R  are  matrices  defined  in  sections  3  and  4  respectively. 
The  matrix  S  is  relevant,  because  by  a  calculation  quite  similar  to 
that  which  led  to  equation  (2.11), 
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^11  0 


Kl3  0 


K22  0 


K 


24 


KT3   0 


K33  K34 


K24   K34   K44 


Kll   0 


Ki3  0 


K22   0 


^24 


^n    0 


K33  0 


^4   0 


^44 


VI 

^    0       ^ 
0 

73 

^ 

.    74     . 

V. 

73+53474 

S  4  37  3+7  4 


Let  X  be  an  eigenvalue  of  S.   Tnen  for  A  /^  1 , 


det 


(l-X)I 


'34 


[   S^3   (l-X)lJ 


det  1 


fl        S34 

[0   (l-X)2l-S,3S34, 


=  det  ((1-X)2l  -  S,.S.,). 


43"^  34' 


We  note  that  X  =  1  is  not  an  eigenvalue  since  83^  and  8^3  are  both 
symmetric,  negative  definite  matrices. 

To  see  that  this  method  is  not  optimal,  we  can  consider  the  case 
when  Ji  is  a  rectangle  cut  in  half.   In  this  case  S^^^  =  S^^^  and  ^^h  = 

\    -...-i'r       ■     ■■-    . 
S43. 

The  eigenvalues  of  S  are  then  of  the  form 

1  +  l/X^d  +  1/2  Rq  +  S^^)). 

By  a  straightforward  calculation  they  are  equal  to  1  +  d^*!'"^  ''.  By 
using  (4.2)  and  (4.3)  one  easily  sees  that  the  eigenvalues  vary  between 
0(l/q  +  1/r)  and  2. 

The  paper  by  Dihn,  Glowinski  and  Periaux  [9],  is  a  report  on  a 
project  to  develop  decomposition  methods  for  nonlinear  problems  in 
fluid  mechanics  e.g.  transonic  flow  on  large  computational  domains. 
Here  we  will  only  discuss  the  two  methods  for  elliptic  problems  which 
were  introduced  in  section  IV  of  their  paper.  In  section  IV. 1,  the 
authors  consider  a  method  in  which  the  unknown  X  is  the  common  value  of 
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3u/9n  on  r  2  .   The  traces  of  the   solutions   on  ^i   and  ^2     ^^^      then 
compared.   The  corresponding  mapping  is,  in  our  notations,  given  by 

and  by  the  results  of  section  5,  this  operator  has  a  condition  number 
which  grows  linearly  with  the  number  of  degrees  of  freedom  associated 
with  Ft  . 

In  section  IV, 2,  the  authors  consider  the  method  which  corresponds 
to  solving  equation  (2.8)  without  preconditioning. 

In  a  recent  paper,  Golub  and  Meyers  [16]  consider  primarily  three 
good  methods  for  the  model  problem  of  section  4  of  this  paper.  Having 
observed  experimentally  that  the  Schur  complement  S  is  quite  close  to  a 
Toeplitz  matrix,  they  derive  an  approximation  of  S  by  solving  a 
simplified  problem  in  which  the  boundaries  of  the  rectangles  fi  ,  and  ^y 
are  moved  to  infinity.  They  also  study  experimentally  the  use  of  the 
preconditioners  J  and  J  (I  +  Jq/4)^'2;  cf.  section  4  of  this  paper. 
They  find  that  the  use  of  the  latter  gives  the  best  rate  of 
convergence. 

The  preconditioners  S^^)  and  J  (I+J^/4)^'2  are  in  fact  quite 
close.  By  slightly  modifying  the  techniques  of  section  4,  we  thus  find 
that 

^%(^  +  J^/4)^/2^  <  x^'^S^^^x  <  (3(r+l))  +  x^^J^CI  +  JpU)^/^x    , 

T 
for  all  XX  =  1. 

In  the  case  considered,  J  and  S^^  commute  and  a  calculation 
shows  that  the  j-th  eigenvalue  of  S^  ■'  exceeds  the  j-th  eigenvalue  of 
J  (I+J^/4)l/2  jjy  ^j^g  factor 

(1  +  (aj^q))2r+2)  /  ^  .  (3(q))2r+2)  , 
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Since  this  factor  approaches  one  very  rapidly,  with  increasing  values 
of  j,  the  matrices  differ  appreciably  only  on  a  low  dimension  subspace. 
The  performance  of  the  two  algorithms  is  in  practice  almost  identical. 

7.    NUMERICAL  EXPERIMENTS  WITH  A  SUBSTRUCTURED  MODEL  PROBLEM 

We  report  on  experiments  with  the  five  point  approximation  of 
Poisson's  equation  on  unions  of  two  rectangles.  The  choice  of  such 
regions  makes  it  feasible  to  conduct  many  experiments  with  very  many 
degrees  of  freedom  since  fast  Poisson  solvers  can  be  used  to  solve  the 
subproblems.  We  note  that  the  five  point  finite  difference 
approximation  results  from  the  use  of  a  conforming  finite  element 
approximation  with  piecewise  linear  basis  functions  on  a  mesh  of  right 
triangles;  see  e.g.  Strang  and  Fix  [32].  Using  the  results  of 
section  4,  we  also  see,  that  on  these  regions,  an  iteration  step  of  our 
algorithms  can  be  carried  out  at  a  speed  which  is  an  order  of  magnitude 
faster  than  a  fast  Poisson  solver  with  general  data.  This  follows  from 
the  fact  that  the  operation  y  =  Sx^  as  well  as  the  preconditioning 
steps,  can  be  carried  out  using  a  few  FFTs.  The  complexity  of  the 
whole  problem  is  therefore  only  about  twice  that  of  the  two  rectangular 
Poisson  subproblems. 

In  the  experiments  reported  here,  we  consider  the  union  of  two 
rectangles  with  corners  at  the  points  (0,0),  (1,0),  (1,1/2),  (0,1/2) 
and  (1/8,1/2),  (k/8,1/2),  (k/8,£/8),  (£/8,Jl/8)  respectively.  We  will 
report  on  cases  when  k  =  3  and  5,  with  Ji  =  6,  8  or  12.  The  mesh  is 
always  uniform  in  both  coordinate  directions  and  the  number  of   points, 

q,  on  r 2  will  therefore  determine  the  discretization  uniquely.   We  have 

2 
used  data  which  are  consistent  with  an  exact   solution  u(x,y)   =  x  + 
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y  -X  e^  cos  y.    We   have   found  no   real  difference  between  the 
performance  of  our  method  for  this  and  other  cases. 

We  first  show,  in  Table  I,  how  five  different  methods  compare  when 
solving  the   problem  with  parameters  k  =   5,  I    =  8   and  using  q  =  63 

unknowns  on  the  interface  ^^  .   We  report  on  the  maximum  error  on  the 

I, 
entire   region  for  each  method  for  the  iteration  at  which  one  of  the 

methods  reaches  the  level  of  the  truncation  error.   The   first  method 

uses   S^^^   as  a  preconditioner ,   the  second  J  while  the  third  method 

considered  solves  a  normal-equation  formulation  of  the  problem,  denoted 

by  Rv[  ,   allowing  the  use  of  R  =  J  as  a  preconditioner.   The  fourth 

column  in  the  table  displays  the  results  when  R  is  used.   This   case 

essentially  shows   the  same  behavior  as  some  of  the  nonoptimal  methods 

analyzed  in  section  6.   The  last  column  shows  the  very  slow  convergence 

when  no  preconditioning  is  used. 


Iteration 

si 

J 

% 

R 

I 

0 

3.73E-1 
1.49E-6 

3.73E-1 
7.82E-5 
1.52E-6 

3.73E-1 
3.63E-3 
2.13E-4 

3.73E-1 
3.95E-2 
1.17E-2 
3.28E-4 

3.73E-1 
1.55E-1 
9.60E-2 
3.78E-2 
1.85E-2 

4 

6 

1.48E-6 

10 

1.48E-6 

1.48E-6 

1.49E-6 

14 

1.48E-6 

1.48E-6 

1.48E-6 

3.17E-6 

Table  1.   Comparison  of  5  different  preconditioners. 
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Table  2  shows  how  the  number  of  Iterations  depends  on  q,  using  the 
two  optimal  preconditioners  S^^-'  and  J.  The  iterations  were  stopped  at 
the  level  of  the  truncation  error  and  the  initial  guess  was  the  zero 
function.  We  note  that  the  over  all  number  of  degrees  of  freedom 
increases  quadratically  with,  q  and  equals  48641  for  q  =  127,  It  is 
clearly  demonstrated  that  both  methods  converge  at  a  rate  which  is 
independent  of  the  mesh  sizeo  The  modest  increase  in  the  number  of 
iterations  reflects  the  increased  accuracy  requirement  as  the 
truncation  error  decreases.  We  employed  the  same  initial  guess  in  all 
these  runs.  _ 


Iterations 

Max.    eirror  oa^ 

q 

s(i) 

J 

s(i) 

_.  ..J 

3 

2 

3 

3.66x10"^ 

3.66x10"^ 

7 

3 

4 

9.59x10-5 

9.55x10-5 

15 

3 

5 

2.45x10-5 

2.43x10-5 

31 

4 

6 

6.09x10"^ 

6.  llx  IQ-^ 

63 

4 

6 

1.49x10"^ 

1.52X  10"^ 

127 

5 

7 

3.02x10-^ 

3.08x10"^ 

Table  2.   Iterations  as  a  function  of  mesh  size. 
The  next  table  shows  a  more  detailed  comparison   between   the   two 
optimal   methods,   giving   the  maximum  error  for  each  method,  at  every 
iteration  for  the  case  when  q  =  127. 
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Iterations 

s(i) 

J 

0 

3.79x10"^ 

3.79x10"^ 

1 

1.25xlO"2 

3.22x10-2 

2 

7.48x10-^' 

4.01x10-3 

• 

3 

2.56x10-5 

,.  5.2^10-^ 

4 

4.42x10"'^ 

8.74x10-5 

5 

3.02x10-7 

1.05x10-5 

1  ,'       -    c           -,■ 

6 

3.02x10^7" 

-1.33x10"^ 

7 

3.02x10-7- 

■"t;  08x1 0-7 

8 

3.02xl0''7-- 

-    3.03x10-7 

Table  3.  A  comparison  of  S*^  -^  and  J. 
Finally,  in  Tables  4  and  5,  we  give  some  information  on  the 
spectra  of  S(S  )-^  and  SJ"  .  We  show  the  two  smallest,  the  fifth  and 
the  two  largest  eigenvalues  for  different  values  of  q  and  for  three 
different  geometries.  The  tables  show  that  the  spectra  of  these 
iteration  operators  are  very  well  behaved  and  are  in  very  close 
agreement  with  the  theoretical  expressions  derived  in  section  4.  We 
note,  in  particular,  that  the  upper  bounds  of  2  and  2/2  appear  in  the 
tables.  The  eigenvalues  in  table  4  have  a  more  pronounced  cluster  and 
they  also  fall  in  a  slightly  smaller  interval  than  those  of  Table  5. 
This  explains  the  faster  convergence  observed  when  S^'-^  is  used.  We 
observe  that  the  eigenvalues  depend  very  weakly  on  the  shape  of  the 
domain,   indicating   that   the   numerical   results  given  in  the  earlier 
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tables  would  not  change  significantly  with   the   domain.   This   is   in 
agreement  with  more  extensive  calculations  that  we  have  performed. 


k   =   5    , 

£    -    6 

k  =   5   , 

Jl    =   8 

k   =   3    , 

1=   12 

X 

q  =  31 

q-63 

q=3i 

q=63 

q  =  31 

q  =  63 

h 

1.714 

1.684 

1.751 

1.713 

1.712 

1.679 

^2 

1.824 

1.776 

1.826 

1.777 

1.820 

1.772 

^5 

1.994 

K985 

1.997 

1.992 

1.996 

1.990 

\-i 

2.000 

2.000 

2.000 

2.000 

2.000 

2.000 

q 

2.000 

2.000 

2.000 

2.000 

2.000 

2.000 

Table  4.   Selected  eigenvalues  of  the  iteration  operator  using  S^^'. 


k   =   5    , 

Jl    =   6 

k   =   5    , 

Jl    =   8 

k   =   3, 

I    =    12 

X 

q  =  31 

q=63 

q  =  31 

q  =  63 

q  =  31 

q  =  63 

h 

1.825 

1.768 

1.778 

1.733 

1.730 

1.692 

X2 

1.868 

1.806 

1.865 

1.804 

1.859 

1.799 

^5 

2.050 

2.014 

2.046 

2.008 

2.046 

2.008 

\-i 

2.822 

2.827 

2.822 

2.827 

2.822 

2.827 

X 

q 

2.827 

2.828 

2.827 

2.828 

2.827 

2.828 

Table  5.  Selected  eigenvalues  of  Che  iteration  operator  using  J. 
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